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j  ABSTRACT 

\J 

Some  methods  are  considered  for  the  estimation  of  probability 
densities.  They  employ  a  linear  approximation  to  either  the  density  or  the 
logistic  density  transform.  The  coefficients  in  the  approximation  are 
estimated  by  maximum  likelihood  and  the  number  of  terms  is  judged  via  an 
information  criterion.  Hence  the  traditionally  difficult  problem  of  judging 
the  degree  of  smoothness  is  carried  out  in  a  relatively  simple  manner. 
Criteria  considered  include  the  penalties  proposed  by  Akaike  and  Schwarz  for 
model  complexity,  together  with  an  empirical  criterion  based  upon  a  plot  of 
the  log-likelihoods.  The  practical  procedures  are  related  to  an  asymptotic 
consistency  argument  and  a  number  of  numerical  examples  are  presented. 
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SIGNIFICANCE  AND  EXPLANATION 


Non-parametric  methods  for  the  estimation  of  probability  densities 
traditionally  experience  severe  technical  and  conceptual  difficulties  in 
judging  the  appropriate  degree  of  smoothness  for  the  estimated  density.  In 
this  paper  some  relatively  simple  procedures  are  discussed  which  employ  either 
polynomial  or  histogram  approximations  and  refer  to  an  information  criterion 
to  judge  the  appropriate  number  of  parameters  in  the  approximation.  The 
duality  between  parameter  parsimony  and  smoothness  is  therefore  exploited. 
Information  criteria  discussed  include  those  due  to  Akaike  and  Schwarz 
together  with  an  empirical  criterion  based  upon  a  plot  of  the  log- 
likelihoods.  The  practical  procedures  are  motivated  by  an  asymptotic 
consistency  argument  and  a  number  of  numerical  examples  are  presented. 
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SOME  PENALIZED  LIKELIHOOD  PROCEDURES  FOR  SMOOTHING 
PROBABILITY  DENSITIES 

Taskin  Atilgan  and  Tom  Leonard 

1.1  Parameter  Parsimony  and  Density  Estimation 

Consider  the  estimation  of  an  unknown  density  f(x),  for  x  6  (a,b) 

from  a  random  sample  of  observations  x1>...,xn.  We  will  employ 

approximations  for  f(x)  taking  the  form 


f  (x)  =  f  (x|9  )  »  A{8  i|>  (x)>  *  X{  J  0  •  (x)} 
m  m  -m  -m-m  .  ,  j  j 

3*1 


(1.1.1 ) 


where  A(»):R+r  is  a  known  transformation,  0  *  (9  ,...,0  )  is  a  vector 

~ni  l  in 

T 

of  unknown  coefficients,  and  ij»  *  (♦,...,♦  )  corresponds  to  a  basis  of 

~m  1  m 

specified  functions,  for  example 


(a)  *  (x)  =  ( 1 ,x,x 


m-1 , 

rX  ) 


(b)  ,  a  suitable  set  of  orthonormal  polynomials 

l  m 

(c)  ♦.(x)  *  I  (x)  for  j  *  1,...,m,  where  the  IT  are  indicator 

J  j  i 

functions  for  a  partition  (Tj,...,T  )  of  the  sample  space.  When 
A(*J  is  the  identity  function  this  leads  to  the  histogram  situation 
discussed  in  section  2.1. 

(d)  corresponding  to  an  appropriate  set  of  B-splines,  so 

i  m 

that  9  ♦  (x)  gives  a  general  formulation  of  a  cardinal  spline  with 
w  m.  in 

m  -  3  knots  (our  ideas  on  splines  will  be  developed  more  fully  in 
another  paper) 

(e)  ♦,...,*  representing  terms  in  a  Fourier  Series 

1  ttl 

The  vector  0  will  be  typically  estimated  by  maximization  of  its  log- 


likelihood 
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n 

L  (0)  -  l  log  A {6  ♦  (x.)}  .  (1.1.2) 

”  "  i-i  '■  -*  1 

A  A 

Let  0  denote  the  maximum  likelihood  vector  and  L  »  L  (0)  the  maximum 
-  mm. 

value  achieved  by  the  likelihood.  These  quantities  will  also  be  important  in 
judging  an  appropriate  value  for  m;  one  possibility  is  to  apply  information 
criteria  proposed  by  Akaike  (1974),  and  Schwarz  (1978),  for  model  selection, 
to  this  situation. 

Akaike  and  Schwarz  suggest  choosing  the  model  which  maximizes  the 
respective  penalized  log-likelihoods 

AIC  -  -  a  (1.1.3) 

and 

SIC  ■  1^  *  |i  log  n  .  (1.1.4) 

These  criteria  introduce  penalties  for  large  values  of  m  and  therefore 
enforce  "parameter  parsimony."  In  the  present  context,  parameter  parsimony 
can  be  viewed  as  representing  smoothness  of  the  unknown  density  f.  The 
penalties  could  therefore  be  referred  to  as  "roughness  penalties”  in  the 
manner  of  I.  J.  Good  (see  section  1.3). 

Akaike  obtained  AIC  by  an  intuitive  argument  based  upon  entropy  measures. 
The  alternative,  SIC,  has  been  frequently  suggested  in  the  Bayesian  literature 
e.g.  by  Jeffreys  (1961),  in  a  variety  of  special  cases.  A  fairly  general 
approximate  Bayesian  formulation  has  been  suggested  by  Leonard  (1981)  and  Chow 
(1981).  The  modification 

SIC'  -  L  -^m  log(n/2* )  (1.1.5) 

might  be  better  for  moderate  sample  sizes. 

The  three  criteria  so  far  described  are  special  cases  of  the  general 
criterion 

GIC  »  L  -  am  (1.1.6) 

in 

where  a  is  the  penalty  per  parameter  in  the  model.  In  this  paper  we  will 
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consider  AIC,  SIC,  SIC*  together  with  an  empirical  procedure  EIC  described 
below  for  choosing  a. 

Note  that  under  fairly  general  regularity  conditions  the  log-likelihood 

Lm  will  possess  the  following  asymptotic  behaviour  as  n  ♦  «* 

-  L  a4.s*  sup  f  f(x)log  f  (x|9  )dx 
n  iq  q  iq  .n 

m 

=  /  f(x)log  f  (x)dx  (1.1.7) 

IQ 

*  n  (f,f  )  , 

Ok  01 

where 

11  (f,f  )  -  -I(f,f  )  +  /  f (x)log  f (x)dx  (1.1.8) 

mm  m 


with 


f  (x)  =  sup  f  (x|8  )  (1.1.9) 

m  m  -m 

where  0  achieves  the  supremum  in  (1.1.7)  and 

«Rt 

I(f,f  )  =  -/  f (x)log{f  (x)/f(x)}dx  (1.1.10) 

IQ  01 

denotes  the  Kullback-Liebler  information  distance  between  f  and  f  .  Note 

m 

that  as  n  ♦  •  maximization  of  n-1Lm  with  respect  to  m  corresponds  to 
minimization  of  this  information  distance. 

Consider,  for  illustrative  purposes,  the  polynomial  situation  described 
above  where  ^(x)  *  (1 ,x,x^, ... ,xm  1 )  and  suppose  that  the  true  density 
takes  the  form 


f(x)  -  fra**x)  *  (1.1.11) 

where  m*  is  the  "true"  number  of  terms  in  the  linear  approximation.  Since 
the  functional  form  of  f  is  often  unspecified  in  practice,  it  will  often  make 
sense  to  simply  refer  to  this  special  form  involving  m*. 

In  this  case,  n(f,fnA)  achieves  the  maximum  over  m  of  the  quantities 
in  (1.1.8),  with  zero  Kullback  information  distance,  yielding  perfectly 
consistent  estimation,  and  we  also  have 


-3- 
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n  (f,f  )  -  n  4(f,f  .)  (1.1.12) 

n  in  di*  mw 

for  m  *  m*,m*  +  1,.... 

since  the  extra  coefficients  for  0  with  m  >  m*  are  set  equal  to  zero  in 

(1.1.9).  Therefore  the  n  will  follow  the  type  of  behaviour  depicted  in 

in 

Figure  1. 


Figure  1 :  Limiting  behaviour  of  log-likelihoods  as  n  *  * 


n 
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Note  that  n  will  tend  to  rise  with  m  until  it  achieves  the  ridge 
m 

at  m  *  m*,  when  the  points  flatten  out.  Therefore  the  problem  at  issue  is 
how  to  infer  this  ridge  point  from  the  observed  quantities  n-1Lm.  We  omit 
detailed  asymptotic  results  which  show  that  as  n  ♦  »,  both  AIC  and  SIC  give 
values  of  m  which  almost  surely  satisfy  m  >  m*,  and  lead  to  consistent 
estimators  for  the  true  density  in  (1.1.10).  Since  SIC  is  based  upon  an 
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asymptotically  larger  penalty  than  AIC  this  suggests  that  SIC  may  lead  to 
values  of  m  which  are  on  the  whole  closer  to  m*  than  those  suggested  by 
AIC. 

The  observed  n-1Lm  will,  for  finite  n,  behave  like  the  theoretical 
points  in  Figure  1,  but  with  random  fluctuation  around  these  points  (see 
Figure  2). 


m 


( -  is  the  upper  convex  boundary) 

When  the  penalty  for  parameter  parsimony  is  a  m,  maximization  of  GIC 
leads  to  the  value  of  m  such  that  the  (nun”1^)  point  in  Figure  2 
maximizes  the  perpendicular  distance  above  the  line  y  ■  a  n  1m.  Therefore 
the  upper  convex  boundary  in  Figure  2  contains  at  its  vertices  all  points 
which  are  optimal  according  to  GIC,  for  some  value  of  a.  For  example,  the 
optimal  AIC  and  SIC  points  must  both  correspond  to  a  vertex  of  this  boundary 


The  upper  convex  boundary  restricts  attention  to  a  convenient  subset  of 
values  of  m.  Sometimes  it  may  be  possible  to  judge  the  ridge-point  m*  in 
Figure  1  by  simple  visual  inspection.  Alternatively  empirical  rules  can  be 
devised  which  estimate  the  position  of  the  ridge  point  e.g.  by  empirically 
choosing  a  value  of  a  rather  than  referring  to  SIC  or  AIC.  An  example  (EIC) 
of  such  a  rule  involves  minimization  with  respect  to  a  and  mQ  of 

e  *  ama  (1.1.13) 

where  ma  >  1  is  a  value  of  m  which  maximizes 

CIO  -  L  -  am  .  (1.1.14) 

SI 

Minimization  of  the  product  in  (1*1.13)  creates  a  compromise  between 
fidelity  to  the  data,  as  represented  by  a  small  value  of  a,  and  smoothness 
of  the  estimate  for  f,  represented  by  a  small  value  for  m.  We  are 
intuitively  speaking  trying  to  obtain  maximum  smoothing  for  minimum  penalty 
am.  Whilst  this  procedure  is  ad  hoc,  our  asymptotic  argument  described  in  the 
next  paragraph  suggests  that  it  gives  a  value  of  m  converging  almost  surely 
to  the  true  value  m*.  Also,  for  finite  n,  it  seems  to  give  reasonable 
answers  in  empirical  situations. 

The  EIC  criterion  gives  an  optimal  value  which  may  be  obtained 
graphically  (see  Figure  2).  We  should  use  the  vertex  of  the  convex  boundary 
which  minimizes  m  times  the  slope  of  the  line  segment  of  the  boundary  lying 
immediately  to  the  right  of  m.  For  the  theoretical  values  in  Figure  1  we  see 
that  this  automatically  leads  to  the  true  value  m*,  where  we  assume  m*  >  1. 
Since  the  n_1Lm  converge  almost  surely  to  the  nm  in  Figure  1,  this 
suggests,  unrigorously,  that  EIC  gives  a  consistent  estimator  for  m*. 
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The  ideas  discussed  above  may  also  be  applied  to  the  histogram  and  spline 
situations  (c)  and  (d)  introduced  at  the  beginning  of  this  section,  but  with 
some  modifications  to  the  limiting  behaviour  of  the  log-likelihoods.  These 
will  be  discussed,  for  histogram  smoothing,  in  section  2.1. 

One  way  of  measuring  the  performances  of  AIC,  SIC,  SIC',  and  EIC  is  via 
simulated  observations  from  a  "true"  density  f.  We  will  judge  the  difference 
between  the  estimated  and  true  values  by  quantities  associated  with  the 
Kullback-Liebler  information  distance  in  (1.1.10). 

Two  contributions  to  the  literature  slightly  related  to  our  approach  are 
by  Crain  (1974),  and  Geisser  and  Eddy  (1979). 

1.2  Further  Uses  of  Parameter  Parsimony 

The  procedures  described  in  section  1.1  depend  upon  the  choice  of  basis 
vector  ♦  (x).  However  they  may  also  be  employed  for  basis  section  e.g.  for 
deciding  whether  to  use,  say,  a  Fourier  Series  approximation  rather  than 
polynomials.  It  is  simply  necessary  to  compare  the  suprema  of  the  penalized 
log-likelihoods  for  the  different  basis  selections.  An  example  will  be 
presented  in  section  2.3. 

All  of  our  criteria  may  be  used  for  estimating  any  one-dimensional 
function  via  a  linear  approximation  e.g.  hazard  functions  for  survival  data, 
the  logit  function  in  non-parametric  bioassay,  and  possible  non-linear 
regression  functions. 

For  the  standard  linear  statistical  model  y.  *  0X.  +  e.  for  j  * 

j  —3  3 

1,...,m,  where  the  e.  are  uncorrelated  normal  errors  with  constant  variance 

3 

2 

a  ,  gic  reduces  to  the  minimization  of 

log(n  R  )  +  am 
z  a 

where  R_  is  the  residual  sum  of  squares.  An  entropy  based  criterion 

B 
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proposed  by  Box  and  Kanamasu  (1973)  uses  a  similar  adjustment,  but  with 
a  *  to  account  for  model  complexity.  It  is  also  possible  to  show  that 
a  *  1  (AIC)  is  equivalent  to  Mallows'  Cp  (see  Mallows,  1973). 

1.3  Other  related  work 

More  formal  approaches  for  the  smooth  estimation  of  curves  have  been 
proposed  by  Wahba  (1977),  Good  and  Gaskins  (1971),  and  Leonard  (1978),  leading 
to  more  sophisticated  roughness  penalties.  For  example,  Wahba  effectively 
modifies  the  log-likelihood  by 

-  -j  A  /  (f(2)(t))2dt  . 

The  parameter  A  controls  the  degree  of  smoothing  and  plays  a  similar 
role  to  m  under  the  present  approach.  This  parameter  is  usually  estimated 
via  cross-validation,  but  the  computations  become  particularly  difficult  in 
non-linear  situations.  An  information  criterion  for  selecting  m  seems  to  be 
much  simpler,  and  seems  to  answer  a  traditinally  difficult  problem  by 
estimating  the  degree  of  smoothness  in  a  relatively  simple  manner. 


2.  APPLICATIONS 


2. 1  Choosing  the  Interval  Width  for  a  Grouped  Histogram 

Let  x1 ,  x2',,,'xn  represent  a  random  sample  of  raw  observations  from  a 
distribution  with  density  f  concentrated,  for  convenience,  on  (0,1).  We 
will  estimate  f  by  a  grouped  histogram  with  m  equal  intervals,  and  will 
judge  the  smoothness  of  the  histogram  by  estimating  m  and  hence  the 
interval  width  m-1 .  This  contrasts  with  a  procedure  proposed  by  Leonard 
(1973)  for  smoothing  the  probabilities  in  a  histogram  with  fixed  intervals. 

Let  n1,...,ng)  denote  the  cell  counts  when  the  n  observations  are  grouped 
into  m  intervals. 

Our  histogram  estimator  has  the  form 

a  m 

f  (x)  -  m  l  P.  I  (x)  (2.1.1.) 

"  j-1  3  Tj 

where  Pj  is  the  proportion  of  raw  observations  lying  in  interval  Tj  = 
(m-1(j-1),  m-1j],  for  j  *  1,...,m,  and  IT  (x)  is  the  appropriate 

j 

indicator  function.  The  log-lifcelihood  for  fixed  m  is  then 


L 

m 


l  log 


i=1 


f  (x.  ) 
m  i 


n 

=  n  log  ra  +  £  n ,  log  P .  . 

j-1  3  3 


(2.1.2) 


Utilizing  a  result  due  to  Basharin  (1959)  we  have,  as  n  ♦  •, 

E(n_1L  )  -  n  (f,f  )  +  —il  +  o^)  (2.1.3) 

n  mm  zt\  d 

n 

where 

E(n  (f,f  )  -  /  f(x)  log  f  (x)dx  (2.1.4) 

n  m  m 

and  f  is  defined  in  a  similar  manner  to  the  quantity  in  (1.1.9). 
m 
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Therefore  (m-1 )/2n  can  be  viewed  as  the  "bias"  when  n_1L.  is  used  to 

in 

estimate  the  expectation  in  (2.1.5).  Note  that  closeness  of  n-1!^  -  (m-1 )/2n 
to  the  true  value  in  (2.1.5)  is  an  indication  that  the  Kullback-Liebler 
information  distance  in  (1.1.10)  is  close  to  zero  i.e.  it  indicates  that  this 
value  of  m  gives  an  estimate  of  f  which  is  close  to  f  according  to  this 
distance  measure. 

Example  1 :  We  generated  n  =  100  observations  from  a  "true"  histogram  with 
m  =  m*  =  5  and  all  probabilities  0.125,  0.25,  0.125,  0.375,  and  0.125. 

With  this  choice  of  f  the  Kullback-Liebler  information  distance  in 
(1.1.10)  is  of  course  minimized  when  m  =  5.  It  yields  a  theoretical  value  of 
0.1153  for  the  expression  in  (2.1.4). 

The  simulation  generated  the  empirical  probabilities  0.11,  0.25,  0.9, 

0.38  when  m  *  5.  Note  from  the  3rd  column  of  Table  1  that  the  bias- 

corrected  log-likelihood  assumes  the  value  11.44,  very  close  to  the  maximum 
theoretical  value  possible  of  11.53.  This  confirms  tht  we  would  like  our 
selection  procedures  to  choose  m  =  5  in  this  situation.  Note,  from  the  3rd 
to  6th  columns,  that  AIC,  SIC,  and  SIC'  are  all  clearly  maximum  when  m  =  5. 

We  have  checked  that  this  is  true  for  all  m  =  1,2,..., 100.  The  three 

criteria  respectively  correspond,  under  GIC,  to  a  =  1,2,3  and  1.38. 

Suppose  in  general  that  the  true  density  f  is  an  equal  interval 
histogram  with  m  =  m*  cells.  Then  the  expression  in  (2.1.4)  will  be 
maximized  when  m  =  m*  and  also  when  m  =  2m*,  3m*,...  .  Therefore,  in 
contrast  to  the  limiting  behaviour  in  Figure  1,  for  the  polynomial  situation, 
the  n  ^Lm,  and  also  the  n  ^1^  -  (m-1 )/2n,  will  possess  the  limiting 
behaviour  depicted  in  Figure  3. 


Figure  3:  Limiting  behaviour  of  histogram  log-likelihoods  as  n 


-*•  <» 


N 


m 


Again  the  empirical  criterion,  EIC,  defined  in  (1.1.13)  and  (1.1.14)  will  lead 
to  a  consistent  estimator  for  ra*  as  n  +  ". 

For  our  simulation  study,  the  log-likelihoods  are  plotted  in  Figure  4  for 
m  »  1,2,..., 20)  they  were  however  calculated  for  all  values  of  m  up  to 
100.  The  vertices  on  the  upper  boundary  are  at  m  =  5,  64,  and  99.  A  visual 
inspection  suggests  that  m  *  5  lies  at  the  beginning  of  the  ridge. 

According  to  GIC,  the  choice  m*  =  5  is  optimal  when  a  lies  between  0.4 
and  0.8.  Higher  choices  are  disregarded  as  too  extreme  to  achieve  optimality 
for  a  reasonable  value  of  a. 

The  primary  choice,  between  m  »  5  and  m  »  64,  can  be  made  by 
reference  to  EIC,  comparing  0.8  *  5  *  4.0  with  0.4  *  64  =  25.6;  again 
suggesting  that  m*  =  5  is  best. 


The  statistician  might  prefer  to  carry  out  a  similar  procedure  but  based 
on  the  plot  of  the  bias-adjusted  log -likelihoods  in  Figure  5.  Note  that 
direct  maximization  amongst  these  adjusted  terms  does  not  in  general  work 
since  it  tends  towards  higher  values  of  m  than  produced  via  an  information 
criterion.  The  latter  again  give  m  *  5  as  best. 

Example  2:  We  next  generated  500  random  observations  from  a  mixture 

0. 3N(-1,1)  +  0.7N(6, 1)  of  two  normal  distributions.  The  results  are 

described  in  Table  2  and  Figure  6.  The  plot  of  bias-adjusted  log-likelihoods 

is  fairly  similar  to  Figure  6  and  is  therefore  omitted.  A  visual  inspection 

of  Figure  6  suggests  that  any  value  of  m  between  11  and  167  would 

probably  be  reasonable.  Note  that  the  quantity  in  (2.1.4)  assumes  the  value 

1013/500  when  f  is  replaced  by  the  true  density  f.  Comparison  with  the 
in 

bias  corrected  log-likelihoods  in  Table  2  suggests  that  m  *  17  is  best,  but 
only  marginally  superior  than  i  *  11, 

In  this  case  both  AIC  and  EIC  gave  m  =  17  and  SIC  and  SIC'  gave  m  = 
11.  In  Pigure  7  the  two  corresponding  histograms  are  compared  with  the  true 
density.  It  seems  that  m  *  17  gives  just  the  right  amount  of  smoothing. 


Table  2:  Values  of  Criteria  for  Example  2 


m 

it  -5=1 

m  2 

AIC 

SIC 

SIC 

1 

-1319 

-1319 

-1320 

-1322 

-1321 

2 

-1276 

-1277 

-1278 

-1282 

-1281 

3 

-1284 

-1285 

-1287 

-1293 

-1291 

4 

-1250 

-1252 

-1254 

-1263 

-1259 

5 

-1137 

-11  39 

-1142 

-1153 

-1148 

6 

-1109 

-1112 

-1115 

-1128 

-1122 

7 

-1090 

-1093 

-1097 

-1112 

-1105 

8 

-1067 

-1071 

-1075 

-1092 

-1085 

9 

-1067 

-1071 

-1076 

-1095 

-1087 

10 

-1048 

-1053 

-1058 

-1079 

-1070 

11 

-1019 

-1024 

-1030 

-1053 

-1043 

12 

-1040 

-1046 

-1052 

-1077 

-1066 

13 

-1026 

-1032 

-1039 

-1067 

-1055 

14 

-1037 

-1044 

-1051 

-1081 

-1068 

15 

-1036 

-1043 

-1051 

-1083 

-1069 

16 

-1020 

-1028 

-1036 

-1070 

-1055 

17 

-1007 

-1015 

-1024 

-1060 

-1044 

18 

-1022 

-1031 

-1040 

-1078 

-1062 

19 

-1022 

-1031 

-1041 

-1081 

-1064 

20 

-1017 

-1027 

-1037 

-1080 

-1061 

21 

-1015 

-1025 

-1036 

-1080 

-1061 

22 

-1007 

-1017 

-1029 

-1075 

-1055 

23 

-1016 

-1027 

-1039 

-1087 

-1066 

24 

-1004 

-1016 

-1028 

-1079 

-1057 

25 

-1009 

-1021 

-1034 

-1087 

-1064 

26 

-1008 

-1020 

-1034 

-1089 

-1065 

27 

-1006 

-1019 

-1033 

-1090 

-1065 

28 

-1008 

-1022 

-1036 

-1095 

-1070 

29 

-1013 

-1027 

-1042 

-1103 

-1076 

30 

-1005 

-1019 

-1035 

-1098 

-1071 

15 
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Figure  7:  Comparison  of  Histograms  and  True  Density 


histogram  obtained 
and  AIC 

histogram  obtained 


by  EIC 
by  SIC,  SIC' 


Example  3:  We  next  consider  the  well-known  chondrite  meteor  data  collected  by 


Ahrens  (1965),  and  previously  analyzed  by  Leonard  (1978),  Silverman  (1981), 
and  Good  and  Gaskins  (1980). 

Table  3:  Silica  Assays  of  22  Chondrite  Meteors 


.04 

.15 

.16 

.78 

.40 

.44 

.45 

.46 

.49 

.54 

.92 

.59 

.64 

.75 

.81 

.33 

.84 

.85 

.87 

.87 

The  log-likelihood  plot  in  Figure  8  yields  the  choice  m  =  16  under  EIC, 
with  a  =  0.85,  a  value  also  selected  by  SIC'  with  a  *  0.63.  However  both 
AIC  and  SIC  choose  m  =  5  with  a  =  1  and  a  =  1.55  respectively.  The 
histograms  for  m  =  5  and  m  =  16  are  depicted  in  Figures  9  and  10. 

Both  histograms  suggest  trimodality  of  the  underlying  distribution, 
agreeing  with  the  results  of  Good  and  Gaskins.  However  m  =  16  seems  to 
create  a  preferable  balance  between  smoothness  and  fidelity  to  the  data. 

For  Examples  1  -  3  we  see  that  EIC  gives  the  most  reasonable  results  in 
all  three  situations  while  AIC,  SIC,  and  SIC'  are  more  variable  in 
performance.  In  the  next  section  EIC  again  fares  well. 


Figure  10:  (m  =  ]6,  SIC' ,  EIC) 


2. 2  Logistic  Density  Estimation 

Apart  from  the  histogram  situation  it  is  not  usually  satisfactory  to  take 

A(»)  in  (1.1.1)  to  be  the  identity  function  since  this  creates  non-negativity 

constraints  on  8  and  it  is  difficult  to  either  estimate  8  by  maximum 
~m  _m 

likelihood  or  to  ensure  that  the  density  integrates  to  unity.  However,  it  is 
sometimes  reasonable  to  follow  Leonard  (1978)  by  considering  the  logistic  form 


~m,m 


.  0  ♦  (s) 

t b  ~m~ra 
fb  e  ds 


(a  <  x  <  b) 


(2.2.1  ) 


as  an  approximation  to  the  underlying  density  f.  For  example,  the  choice 

♦  (x)  *  (x,x2, . . . ,xm)T  ensures  that  the  first  m  sample  moments  constitute  a 
in 

set  of  sufficient  statistics  for  8  .  Transformation  to  a  set  of  suitable 

win 

orthogonal  polynomials  yields  well-conditioned  Hessians  for  the  iterative 
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solution  of  the  likelihood  equations  by  Newton-Raphson.  The  latter  involves 


maximization  of  the  log- likelihood  function 


where 


L  {6  )  =  9  t  -  n  log  exp{0  ^  (x)}dx 
m  ,m  -m-m  b  -m-m 


t  -  l  *<x.  ) 
m  ,=1  -n  i 


(2.2.3) 


(2.2.4) 


and  can  be  tackled  using  standard  optimization  packages. 

In  our  numerical  example  we  completed  this  maximization  for  m  * 

1,2/.. .,10  (by  which  time  we  had  obviously  reached  the  ridge-point)  and  then 
selected  m  via  our  information  criteria. 

Figure  11:  Log-Likelihoods  for  Logistic  Example 
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Example:  200  random  numbers  are  generated  from  the  mixture  0.3N(-1,1)  + 
0.7N( 6, 1 )  of  two  normal  distributions.  Using  a  basis  of  polynomials  up  to 
degree  m  we  calculated  the  log-likelihoods  depicted  in  Figure  11. 

A  visual  inspection  suggests  that  either  ra  *  4  or  m  *  6  is 
appropriate.  However  EIC,  SIC,  SIC',  and  AIC  respectively  chose  m  *  4,  6, 
10,  and  10,  with  a  -  3.21,  2.5,  1.7,  and  1.0.  Figures  12  -  14  compare 
the  corresponding  density  curves  with  the  true  curve  and  a  histogram  of  the 
raw  observations.  It  seems  that  EIC  and  m  *  4  gives  about  the  right  amount 
of  smoothing  for  a  small  number  of  terms  in  the  approximation. 
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2.3  Basis  Selection 


It  is  possible  to  use  our  penalized  log-likelihoods  to  compare  different 
choices  of  basis  functions.  In  the  following  example  we  illustrate  this 
approach,  as  part  of  a  more  general  simulation  study,  by  comparing  Hermite  and 
Laguerre  functions. 

We  completed  500  repeated  simulations,  and  each  time  simulated  a  random 
sample  of  size  100  from  a  chi-squared  distribution  with  14  degrees  of 
freedom.  In  each  case  we  used  the  logistic  method  described  in  section  2.2, 
to  estimate  the  density,  and  the  500  repetitions  enabled  us  to  get  an  accurate 
estimate  (IMSE)  together  with  its  standard  deviation  (STND),  and  average 
values  for  the  log- likelihoods  together  with  AIC,  SIC,  and  SIC1.  The  results 
are  summarized  in  Tables  4  and  5. 

From  the  values  of  m  which  we  were  able  to  calculate  we  see  that  the 
preference  for  Laguerre  functions  is  very  clear  under  IMSE,  AIC,  SIC,  and 
SIC',  giving  useful  insight  which  would  be  difficult  to  judge  intuitively. 
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Table  4: 

Simulated 

Results 

(Hermite 

Functions ) 

m 

IMSE 

STND 

Lm 

AIC 

SIC 

SIC' 

2. 

.01206 

.00066 

-3113 

-3115 

-3119 

-3117 

4. 

.01206 

.00066 

-2819 

-2823 

-2831 

-2828 

6. 

.01207 

,00066 

-2560 

-2566 

-2578 

-2573 

8. 

.01207 

.00066 

-2341 

-2349 

-2366 

-2359 

10. 

.01208 

.00066 

-2116 

-2126 

-2147 

-2138 

12. 

.01209 

.00065 

-1964 

-1976 

-2002 

-1991 

14. 

.01209 

.00065 

-1820 

-1834 

-1863 

-1850 

Table  5:  Simulated  Results  (Laguerre  Functions) 


m 

IMSE 

STND 

Si 

AIC 

SIC 

SIC 

2. 

.01127 

.00071 

-584 

-586 

-538 

-587 

4. 

.00391 

.00042 

-301 

-305 

-310 

-307 

6. 

.00220 

.00023 

-276 

-282 

-290 

-285 

8. 

.00212 

.00025 

-273 

-281 

-291 

-284 

10. 

.00211 

.00026 

-270 

-280 

-293 

-284 

12. 

.0021 1 

.00026 

-276 

-288 

-303 

-292 

14. 

,00211 

,00026 

-264 

-278 

-296 

-283 
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